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Rough-wall  channel  analysis  using  suboptimal 

control  theory 

By  O.  Flores  I,  J.  Jimenez  X and  J.  Templeton 


1.  Introduction 

The  original  aim  of  this  work  was  to  shed  some  light  on  the  physics  of  turbulence 
over  rough  walls  using  large-eddy  simulations  and  the  suboptimal-control  wall  boundary 
conditions  introduced  by  Nicoud  et  al.  (2001).  It  was  hoped  that,  if  that  algorithm 
was  used  to  fit  the  mean  velocity  profile  of  the  simulations  to  that  of  a rough-walled 
channel,  instead  of  to  a smooth  one,  the  wall  stresses  introduced  by  the  control  algorithm 
would  give  some  indication  of  what  aspects  of  rough  walls  are  most  responsible  for  the 
modification  of  the  flow  in  real  turbulence.  It  was  similarly  expected  that  the  structure 
of  the  resulting  velocity  fluctuations  would  share  some  of  the  characteristics  of  rough- 
walled  flows,  thus  again  suggesting  what  is  intrinsic  and  what  is  accidental  in  the  effect 
of  geometric  wall  roughness. 

A secondary  goal  was  to  study  the  effect  of  ‘unphysical’  boundary  conditions  on  the 
outside  flow  by  observing  how  a relatively  major  change  of  the  target  velocity  profile, 
and  therefore  presumably  of  the  applied  wall  stresses,  modifies  properties  such  as  the 
dominant  length  scales  of  the  velocity  fluctuations  away  from  the  wall. 

As  will  be  seen  below,  this  secondary  goal  grew  more  important  during  the  course  of 
the  study,  which  was  carried  out  during  a short  summer  visit  of  the  first  two  authors 
to  the  CTR.  It  became  clear  that  there  are  open  questions  about  the  way  in  which 
the  control  algorithm  models  the  boundary  conditions,  even  for  smooth  walls,  and  that 
these  questions  make  the  physical  interpretation  of  the  results  difficult.  Considerable 
more  work  in  that  area  seems  to  be  needed  before  even  relatively  advanced  large-eddy 
simulations,  such  as  these,  can  be  used  to  draw  conclusions  about  the  physics  of  wall- 
bounded  turbulent  flows. 

The  numerical  method  is  the  same  as  in  Nicoud  et  al.  (2001).  The  modifications  in- 
troduced in  the  original  code  are  briefly  described  in  §2,  but  the  original  paper  should 
be  consulted  for  a full  description  of  the  algorithm.  The  results  are  presented  in  §3  and 
summarized  in  §4.  The  elementary  properties  of  turbulence  over  rough  walls  which  are 
used  in  the  text  have  been  taken  from  recent  reviews  such  as  Raupach  et  al.  (1991)  or 
Jimenez  (to  appear  in  2004). 


2.  Simulations 

Turbulent  channel  flow  has  been  simulated  using  the  code  presented  by  Nicoud  et  al. 
(2001)  with  ReT  = urh/v  = 1000,  where  uT  is  the  friction  velocity  and  h is  the  channel 
half-width.  The  coordinates  used  are  x,  y and  z,  indicating  streamwise,  wall-normal  and 
spanwise  directions.  The  streamwise  and  spanwise  periodicity  lengths  are  Lx  = 2i r and 
Lz  = 2tt/3  for  all  the  cases  presented  here.  Except  when  explicitly  noted  to  the  contrary, 

f School  of  Aeronautics,  Universidad  Politecnica,  28040  Madrid,  Spain 
1 Also  School  of  Aeronautics,  Universidad  Politecnica,  28040  Madrid,  Spain 
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Rer 

Ai+ 

A y+ 

A z+ 

Lx/h 

Lz/h 

Nx  x Nv  x Nz 

a+ 

A U+ 

SI 

1000 

196 

61 

65 

2n 

2tt/3 

32  x 33  x 32 

6 x 10-2 

0 

S2 

1000 

196 

61 

65 

27T 

2tt/3 

32  x 33  x 32 

1.5  x 10"2 

0 

R1 

1000 

196 

61 

65 

2n 

2tt/3 

32  x 33  x 32 

6 x 10"2 

9 

Table  1.  Test  cases.  The  subindices  x , y and  z indicate  the  streamwise,  wall-normal  and  span- 
wise  directions.  The  N’s  are  the  number  of  grid  points  in  each  direction,  the  A’s  are  the  cor- 
responding grid  spaces  and  the  L' s are  box  lengths.  AU+  is  the  roughness  function,  defined  in 
(3.7),  and  a+  is  the  parameter  in  equation  (2.1). 


all  the  variables  are  normalized  with  uT  and  h.  The  superindex  + is  used  for  variables 
expressed  in  wall  units. 

The  characteristics  of  the  three  main  simulations  used  in  this  paper  are  summarized 
in  table  1.  The  letter  S indicates  smooth  channels,  and  R is  used  for  rough  ones.  The 
reference  velocity  profile  used  for  SI  and  S2  is  obtained  from  a direct  simulation  by  del 
Alamo  et  al.  (2003)  with  ReT  = 950.  The  rough-wall  profile  is  described  in  the  next 
section. 

The  LES  code  uses  a standard  dynamic  Smagorinsky  subgrid-scale  stress  model  with- 
out explicit  grid  filtering.  The  spatial  discretization  is  a second-order  finite-difference 
scheme  on  a staggered  grid,  and  the  time  integration  is  third-order  Runge-Kutta  with 
an  implicit  scheme  for  the  wall-normal  viscous  terms.  The  flow  is  driven  by  a constant 
mean  pressure  gradient  that  balances  the  wall  stresses,  dxp  + — — 1.  The  overbax  stands 
for  averaging  over  wall-parallel  planes,  while  { ) will  be  reserved  for  temporal  averaging. 
With  the  single  exception  described  at  the  end  of  §3.1,  the  time  step  At  is  constant. 

The  boundary  conditions  in  the  streamwise  and  spanwise  directions  are  periodic,  and 
the  impermeability  condition  is  imposed  at  the  walls.  The  two  other  boundary  conditions 
are  the  total  shear  stresses  rxy  and  rzy  at  each  point  of  each  wall.  They  are  adjusted  at 
each  time  step  to  minimize  a cost  function  J,  which  measures  the  difference  between  the 
plane- aver  aged  velocity  and  a given  mean  velocity  profile  Uref. 

= J [(U  - Uref)2  + W2}  dy  + Q + 0g)  , (2.1) 

where  the  control  variables  are  vectors  of  size  2 x Nx  x Nz 

<t>u  = {rxy  \y=2 h,  -Txy  |„=o),  <t>w  = (t“  |y=2fc,  ~Tzy  ll/=o),  (2.2) 

and  a is  a parameter  with  dimensions  [a]  = L/U 2 that  will  be  discussed  in  the  next 
section.  Primed  variables  refer  to  fluctuations  with  respect  to  the  plane  average,  so  that 
<t>'  = <j>  — cj>.  Note  that  the  global  conservation  of  momentum  ensures  that  </>  satisfies 

(Wy)  \y~2h  -(7%)  |y=0=  1,  (2.3) 

and 

(7%)  ly=2h  ~(7fy)  ly=0—  0,  (2.4) 

but  that  the  instantaneous  mean  stresses  at  each  wall  may  oscillate  in  time. 

The  minimization  of  J implies  an  iterative  scheme  in  which  an  adjoint  problem  is 
solved  in  order  to  compute  an  estimation  for  the  gradient  of  the  cost  function.  The 
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Figure  1.  Convergence  history  of  the  mean  velocity  at  the  first  plane  off  the  wall,  for  a 
smooth-walled  channel. , Using  the  correction  (2.5); , without  the  correction. 


algorithm  is  suboptimal  in  the  sense  that,  instead  of  minimizing  the  averaged  value  of  J 
over  a long  period  of  time,  the  optimization  is  computed  over  each  individual  time  step. 

Several  modifications  were  introduced  in  the  code  during  preliminary  attempts  to 
match  smooth-wall  velocity  profiles  at  different  Reynolds  numbers.  The  most  important 
was  probably  the  use  of  the  more  efficient  Brent’s  method  for  the  optimization  scheme 
(Press  et  al.  1993),  instead  of  the  simpler  relaxation  used  by  Nicoud  et  al.  (2001).  This 
allowed  the  stable  utilization  of  the  algorithm  over  a wider  range  of  parameters. 

Another  modification  involves  the  recalculation  of  the  mean  wall  stress.  In  the  original 
code,  the  value  of  + = u+2  was  adjusted  at  each  time  step  so  that  the  mean  velocity 
in  the  first  plane  off  the  wall  satisfied  the  desired  logarithmic  law, 

-gi-ilog(j/t=1u+)-A  = 0,  (2.5) 

The  idea  was  that  the  mean  stress  is  directly  determined  by  the  averaged  momentum 
balance,  and  does  not  have  to  be  controlled.  In  each  time  step  the  optimum  wall  stresses 
were  computed  using  the  control  algorithm,  and  their  mean  value  was  later  corrected  to 
the  u2  obtained  from  (2.5).  Notice  that  uf  = 1 satisfies  (2.5)  whenever  ut=1  is  given  by 
the  desired  logarithmic  law.  While  the  general  argument  is  sound  in  a timed-averaged 
sense,  this  procedure  ‘second-guesses’  the  control  algorithm,  since  the  integrated  mo- 
mentum equation  is  already  incorporated  in  the  code,  and  its  effect  was  found  to  be 
pernicious.  The  convergence  of  the  original  code  was  oscillatory  and  the  time  history  of 
most  variables  suggested  a poorly-damped  instability  of  the  control  algorithm.  Removing 
the  correction  step  suppressed  the  oscillations  and  the  instability.  Two  typical  conver- 
gence histories,  with  and  without  the  stress  correction,  are  shown  in  figure  1.  All  the 
results  below  are  obtained  without  using  (2.5). 


3.  Results 

3.1.  The  cost  function 

A point  that  needs  some  discussion  is  the  form  of  the  cost  function  (2.1).  The  second 
term  on  the  right-hand  side  of  (2.1)  represents  the  energy  of  the  control  variables  4>u  and 
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FIGURE  2.  Profiles  of  (a)  the  mean  velocities  and  (6)  the  streamwise  velocity  fluctuations  for  a 
smooth-walled  channel  at  ReT  ss  1000,  using  different  cost  functions,  o , case  SI;  • , case  S2; 
, direct  simulation  (del  Alamo  et  al. 2003). 

4>w , and  is  introduced  to  ensure  numerical  stability.  The  proportionality  coefficient  a is 
essentially  arbitrary  from  a physical  point  of  view. 

The  value  of  a used  here  for  the  cases  SI  and  R1  is  the  one  given  by  Nicoud  et  al. 
(2001)  as  being  the  smallest  one  that  kept  their  simulation  stable.  The  effect  of  a larger  a 
should  be  to  lower  the  variance  of  tw,  weakening  the  strength  of  the  control,  and  therefore 
presumably  degrading  the  quality  of  the  resulting  mean  profile.  It  is  nevertheless  possible 
to  ‘overcontrol’  a system,  in  the  sense  that  some  properties  which  are  not  taken  into 
account  in  the  cost  function  may  be  contaminated  by  the  control  itself.  It  was  noted  by 
Nicoud  et  al.  (2001)  that  the  velocity  fluctuations  in  the  simulated  channels  are  stronger 
than  those  in  direct  simulations  or  in  experiments,  especially  near  the  wall.  Since  those 
fluctuations  are  driven  by  the  wall  stresses,  which  are  in  turn  influenced  by  a,  it  is 
conceivable  that  the  reason  for  the  high  fluctuations  is  that  a was  chosen  too  small. 

A measure  of  the  strength  of  the  control  is  the  temporal  variability  of  u,  compared 
with  its  ‘expected’  variation.  If  the  velocities  are  considered  to  be  uncorrelated  random 
variables,  which  is  probably  not  too  far  from  the  truth  for  a coarse  LES,  the  expected 
variance  of  u would  be 

4 = ( NXNZ )-1^,  (3.1) 

which  has  to  be  compared  with  the  actual  measured  value 

a2  = ((u-(u))2).  (3.2) 

The  ratio  o /a p should  be  of  order  unity  in  a natural  flow,  and  its  deviation  from  unity 
is  a measure  of  the  strength  of  the  control. 

It  is  not  clear  which  is  the  right  normalization  for  a,  but  the  observation  in  Nicoud  et  al. 
(2001)  that  the  effect  of  the  control  is  mainly  felt  on  the  velocities  near  the  wall  suggests 
that  h is  not  a relevant  length  scale,  and  that  a more  natural  normalization  would  be 
wall  units.  The  value  recommended  by  Nicoud  et  al.  (2001)  is  then  a+w6x  10-2,  and 
results  in  a /crp  «2x  10-3.  Both  the  small  magnitude  of  the  dimensionless  a and  of  the 
ratio  of  the  standard  deviations  strongly  suggest  that  the  system  is  overcontrolled. 

A detailed  analysis  of  the  effect  of  a is  beyond  the  scope  of  this  work,  and  it  was 
not  attempted,  but  a simple  approximation  of  the  cost  function  can  be  used  to  estimate 
the  relevant  trends.  Assume  that  the  velocity  u near  the  wall  is  essentially  a (random) 
function  of  the  local  wall  stresses,  u — f(rw),  where  all  the  equations  in  the  following 
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argument  have  to  be  understood  as  being  applied  at  a given  x and  2.  We  can  then 
approximate  the  cost  function  as, 

3 « 2 8(u  - Uref)2  + 2 ttpf,  (3.3) 


where  S is  a measure  of  the  height  over  which  the  velocity  is  modified  by  the  control, 
and  u is  taken  at  some  representative  distance  from  the  wall  of  order  <5.  The  factor  of 
two  reflects  the  contribution  from  both  walls.  The  minimization  of  J requires  that  its 
derivative  with  respect  to  all  the  tw  should  vanish,  giving  for  each  grid  point 


dr ” NXNZ  drw  NXNZ 


= 0, 


(3.4) 


where  we  have  assumed  that  t“  does  not  depend  of  the  control.  Assuming  that  f(rw) 
can  be  linearized  around  r57  and  expressing  u in  terms  of  f(rw),  we  obtain  after  some 
algebra  that 

T”'  = {Uref-J^}  6f^ (3.5) 

1 + 8f? /a 

where  fr  is  d//dr“  evaluated  at  r“.  This  is  an  interesting  expression  because  it  captures 
the  tendency  of  t'w  to  increase  with  decreasing  a,  but  saturates  to  some  finite  value  when 
a — > 0.  In  that  limit  the  cost  function  becomes  u — Uref  oc  a. 

This  has  the  simple  interpretation  that,  if  there  is  a solution  to  the  unconstrained 
control  problem,  i.e.  if  there  is  a tw  such  that  the  first  term  in  (2.1)  vanishes,  the  control 
algorithm  approximately  identifies  it  for  some  small  value  of  a,  and  any  further  decrease 
in  a does  not  change  that  solution  appreciably.  If  such  solution  does  not  exist,  t'w  and 
u'  would  diverge  as  when  a — > 0,  and  the  cost  function  would  saturate  to  some 
non-zero  minimum  value.  This  would  appear  in  our  simple  model  as  ft  — 0. 

We  tested  this  behavior  by  running  case  S2,  in  which  a was  divided  by  four  with 
respect  to  the  two  other  cases.  The  new  value  of  a decreases  cr/crp,  as  expected,  while 
the  mean  velocity  profile  improves  everywhere  (figure  2 a),  and  the  intensity  of  the  velocity 
fluctuations  near  the  wall  increases  (figure  26).  All  this  is  in  qualitative  agreement  with 
the  previous  analysis.  The  ratio  between  the  standard  deviations  of  the  wall  stresses  in 
cases  S2  and  SI  is  2:1,  which  is  smaller  than  the  ratio  of  4:1  which  would  correspond  to 
an  a-1  behavior,  and  suggests  some  degree  of  saturation  towards  an  exact  solution. 

To  test  whether  such  a solution  exists,  a second  experiment  was  run  with  a = 0.  The 
velocity  fluctuations  grew  still  further,  apparently  without  bound,  leading  to  numerical 
instability  for  any  given  constant  time  step.  When  the  code  was  modified  to  run  at  a 
constant  CFL,  we  were  able  to  trace  the  growth  of  the  fluctuating  velocities  in  the  first 
plane  off  the  wall  to  values  of  the  order  of  u'+  ~ 600,  implying  that  the  suboptimal 
boundary  conditions,  at  least  in  their  present  form,  are  not  able  to  drive  the  LES  to  full 
agreement  with  the  desired  profile.  The  errors  in  figure  2(a)  suggest  that  the  difficulty  is 
not  with  the  behavior  of  the  near-wall  region,  but  with  the  center  of  the  channel.  It  is  an 
interesting,  although  unresolved,  question  whether  this  is  a limitation  of  the  LES  model 
itself,  which  cannot  reproduce  the  phenomena  responsible  for  the  wake  component,  or  of 
the  boundary  conditions. 


3.2.  Rough  walls 

The  most  important  effect  of  roughness  on  the  mean  velocity  above  the  buffer  zone  is 
a constant  velocity  decrement,  A U+  (Raupach  et  al.  1991;  Jimenez,  to  appear  2004). 
A second  effect  is  a shift  A y in  the  wall-normal  coordinate,  due  to  the  uncertainty  in 
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Figure  3.  (a)  Mean  velocity  profiles.  (6)  Streamwise  velocity  fluctuations,  (c)  Wall-normal 

velocity  fluctuations,  (d)  Spanwise  velocity  fluctuations. , smooth  wall  DNS,  ReT  = 950; 

, reference  velocity  profile  for  the  rough  case;  o , case  SI;  a , case  Rl. 


the  position  of  the  wall.  In  our  case,  in  which  the  velocity  is  only  computed  at  fairly 
large  values  of  y+ , the  wall- normal  shift  is  negligible  and  the  mean  velocity  profile  can 
be  written  as 

U+(y+)  = - log(t/+)  + 8.5  - - log(fc+ ) + -n{y/h),  (3.6) 

K /C  K 

or 

U+(y+)  = - log(y+)  + 5.1  - AU+  + -Il(y/h),  (3.7) 

were  II  is  the  wake  function,  and  either  A U+  or  kf  characterize  the  effect  of  roughness 
on  the  mean  velocity  profile. 

Using  equation  (3.7)  we  construct  our  rough  mean  velocity  profile  for  case  Rl  just 
subtracting  a constant  AU+  = 9+  from  the  smooth  mean  velocity  profiles  in  del  Alamo 
et  al.  (2003).  This  corresponds  to  a fully-rough  flow  with  an  equivalent  sand  roughness 
kf  « 140. 


3.3.  Statistics 

Figure  3 shows  the  mean  velocity  and  the  velocity  fluctuations  profiles  for  cases  SI  and 
Rl.  The  mean  velocities  share  some  characteristics  with  those  obtained  by  Nicoud  et 
al.  (2001),  even  though  the  latter  used  as  reference  profile  a logarithmic  law  without  a 
wake  function.  The  error  between  the  mean  velocity  and  Uref  is  given  in  figure  4 for 
the  three  cases  computed.  It  is  quite  small  near  the  wall,  but  grows  to  0(1)  at  the 
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Figure  4.  Error  in  the  mean  velocity  profiles  with  respect  to  their  references,  o , case  SI;  • 

case  S2;  a , case  Rl. 


center  of  the  channel.  Between  the  second  and  the  third  grid  point  there  is  a weak 
discontinuity  that  changes  the  intercept  of  the  logarithmic  law.  Its  origin  is  probably  the 
spatial  discretization,  because  it  appears  at  the  same  grid  point  in  Nicoud  et  al.  (2001) 
even  if  their  resolution  is  coarser  than  ours  by  a factor  of  four. 

The  root-mean-squared  velocity  fluctuations  are  presented  in  figures  3(6),  3(c)  and 
3(d).  The  streamwise  and  spanwise  fluctuations  near  the  wall  are  higher  in  all  the  large- 
eddy  simulations  than  in  the  direct  simulation,  although  the  agreement  is  quite  good  in 
the  core  region  for  the  streamwise  direction.  The  wall-normal  and  spanwise  fluctuations 
are  slightly  underpredicted  by  the  suboptimal  code  in  the  core  region.  The  high  values 
for  the  streamwise  velocity  fluctuations  near  the  wall  and  the  velocity  fluctuations  in 
the  core  region  agree  with  the  results  of  Nicoud  et  al.  (2001),  but  the  exceptionally  high 
values  of  the  spanwise  velocity  fluctuations  near  the  wall  shown  in  figure  3(d)  do  not. 

The  lower  values  of  the  velocity  fluctuations  near  the  wall  in  case  Rl  with  respect  to 
SI  are  consistent  with  the  effect  of  roughness  on  physical  turbulence.  This  is  often  inter- 
preted as  the  interference  of  the  roughness  elements  with  the  smooth-wall  self-sustaining 
turbulence  cycle  (Jimenez  & Moin  1991),  as  explained  by  Jimenez  (to  appear  in  2004). 
This  interpretation  is  unlikely  in  the  present  case  because  the  resolution  is  too  coarse 
in  all  three  directions  to  capture  the  scales  associated  with  the  regeneration  cycle,  and 
the  most  likely  explanation  is  that  the  decrease  of  the  fluctuations  is  due  to  the  lower 
turbulence  production  near  the  wall  due  to  the  weaker  velocity  gradients  across  the  first 
grid  element  of  the  rough  reference  profile. 

The  wall  stresses  T™y  and  r™y  provided  by  the  control  algorithm  are  studied  using  their 
probability  density  functions  (p.d.f.),  which  are  shown  in  figure  5.  They  are  compared 
at  y+  = 60  with  the  p.d.f. ’s  of  the  shear  stresses  in  a DNS  channel  with  ReT  = 550  (del 
Alamo  et  al.  2003),  after  box-averaging  the  latter  to  the  LES  grid  (A*  x Az  = 196+x65+). 
Although  it  was  not  possible  to  post-process  the  shear  stresses  from  the  ReT  = 950  in  time 
for  this  report,  there  is  probably  little  differences  between  them  an  those  at  Rer  = 550, 
because  this  part  of  the  flow  scales  approximately  in  wall  units. 

It  is  clear  from  figure  5(a)  that  the  r™y  in  the  suboptimal  simulations  axe  different  from 
the  DNS  one,  especially  in  that  the  former  are  not  able  to  reproduce  the  asymmetry  of 
the  p.d.f.  and  try  to  reproduce  the  mean  value  of  the  shear  stresses  by  displacing  their 
modes  to  positive  values.  As  expected,  the  standard  deviations  behave  in  the  same  way 
as  the  velocity  fluctuations,  Rl  < SI  < S2.  Figure  5(c)  and  5(d)  show  the  p.d.f. ’s  of  the 
three  cases  when  the  stresses  are  normalized  with  theirs  mean  value  (r)  and  standard 


420  0.  Flores,  J.  Jimenez  & J.  Templeton 


FIGURE  5.  Probability  density  functions  (p.d.f.)  of  the  shear  stresses  at  the  wall.  The  shear 
stresses  in  figures  (a)  and  ( b ) are  expressed  in  wall  units,  while  in  figures  (c)  and  (d)  they  are 

normalized  as  in  equation  (3.8). , DNS  (see  text  for  details); , Gaussian;  o , case  SI; 

• , case  S2;  A , case  Rl. 


deviation  aT, 

rN  r-  (r) 
ar 

The  collapse  of  the  three  curves  on  the  Gaussian  distribution,  and  the  differences  with 
the  DNS  results,  suggest  that  there  is  little  physical  information  on  the  shear  stresses 
given  by  the  control  algorithm. 

Another  way  of  characterizing  the  shear  stresses  given  by  the  control  is  the  study  of 
their  spectral  distribution.  We  can  define  the  spectrum  Ey 

Eij(kx,kz)  = (Tij?*j),  (3.9) 

where  t, ij(kx,  kz)  are  the  Fourier  coefficients  of  the  two-dimensional  Fourier  transform 
of  the  wall  stresses  T^j.  The  asterisk  * indicates  complex  conjugation  and  kx,kz  are  the 
wave  numbers. 

Figure  6 shows  these  premultiplied  spectra  for  the  cases  SI  and  Rl,  as  functions  of 
the  wavelengths  \x  = 2ir/kx  and  Xz  — 2n/kz.  Because  each  spectrum  is  normalized  with 
its  maximum,  it  only  contains  information  about  the  wavelengths.  There  is  very  little 
differences  between  the  two  cases. 

It  is  important  to  notice  that  the  most  energetic  modes  are  located  in  the  large- 
wavelength  end  of  the  spectra,  specially  as  regards  their  widths.  The  ‘infinitely  wide’ 
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FIGURE  6.  Two-dimensional  premultiplied  spectra  of  the  stresses  at  the  wall  as  functions  of  the 
streamwise  and  spanwise  wavelengths,  (a)  Exy,  ( b ) Ezy.  Shaded  contours  are  from  the  SI  case 
and  lines  are  from  Rl.  Each  spectrum  is  normalized  with  its  maximum. 


(a,  0)  modes  for  the  r™y  spectrum  sums  up  to  26%  of  the  energy  in  the  smooth  case  and 
to  35%  in  the  rough  one,  but  the  energy  in  the  ‘infinitely  long’  (0,  /?)  modes  is  below  1% 
of  the  total  in  both  cases.  On  the  other  hand,  for  the  r™y  spectrum  the  energy  contained 
in  the  (a,  0)  modes  is  negligible,  while  the  energy  contained  in  the  (0,  (3)  modes  sums  7% 
of  the  total  energy  in  the  smooth  case  and  9%  in  the  rough  case.  That  suggests  that  the 
spectrum  of  r™y  is  wider,  but  not  too  much  longer  than  the  simulation  box,  while  the  r™y 
spectrum  almost  fits  in  it.  The  behavior  of  this  quantity  in  the  DNS  is  not  known. 

It  is  interesting  to  analyze  how  the  control  variable  characteristics  affect  flow  variables 
such  as  the  velocities.  In  order  to  do  that,  we  have  studied  the  energy  spectrum  of  the 
velocity  components,  Euu,  Evv  and  Eww,  defined  as 

Euu{hx,  kz)  = {'akx,kzukx,kz'}T  (3.10) 

The  corresponding  one-dimensional  spectra  axe  obtained  adding  equation  (3.10)  over  a 
certain  direction  in  Fourier  space. 

Figures  7(a)  and  7(b)  show  the  one-dimensional  premultiplied  spectrum  of  the  stream- 
wise  velocity  components  at  y+  — 30,  which  is  the  first  grid  point  in  the  mesh.  Although 
the  gradient  of  this  variable  is  determined  directly  by  the  control,  the  agreement  between 
the  suboptimal  cases  and  the  DNS  results  is  reasonable  good,  at  least  as  much  as  can  be 
expected  from  the  coarse  grid  being  used.  The  same  happens  for  the  wall-normal  velocity 
component  (not  shown  here),  but  not  for  the  spanwise  component,  whose  spectrum  is 
shown  in  figures  7(c)  and  7(d).  The  spanwise  velocity  has  a very  energetic  mode  which 
spans  the  full  width  of  the  box  and  almost  all  of  its  length.  When  the  instantaneous 
velocity  fields  are  examined,  we  observe  fairly  regular  diagonal  bands  of  positive  or  neg- 
ative spanwise  velocity.  This  is  reminiscent  of  the  instabilities  observed  by  Jimenez  et  al. 
(2001)  in  a channel  with  a porous  wall,  where  they  were  traced  to  a coupling  between 
a weakly-damped  mode  of  the  impermeable  channel  and  the  porous-wall  condition.  In 
that  case  the  result  was  a series  of  spanwise  rollers  spanning  the  full  height  and  width 
of  the  channel,  and  it  is  conceivable  that  a similar  coupling  might  result  in  the  present 
structures. 

The  spurious  perturbation  of  the  spanwise  velocity  disappears  as  we  move  a few  grid 
points  away  from  the  wall,  as  can  be  observed  in  the  spectra  in  figures  8.  This  relatively 
fast  relaxation  away  from  the  wall  of  the  defects  of  the  wall  layer,  which  recalls  the  similar 
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Figure  7.  One-dimensional  premultiplied  velocities  spectra  at  y+  — 30,  normalized  in  wall 
units  and  plotted  as  functions  of  the  wavelengths  Xx  (a,  c)  and  Az  (b,  d).  (a)  and  (b),  streamwise 
component;  (c)  and  (d)  spanwise  component. , DNS;  o , case  SI;  • , case  S2;  a , case  Rl. 


relaxation  of  the  fluctuation  intensities  in  figures  2(b)  and  3(b),  suggests  that  the  effect 
of  the  wall  is  relatively  local  to  small  values  of  y , and  that  the  main  role  of  the  boundary 
conditions  is  to  provide  a correct  intercept  for  the  mean  velocity. 


4.  Conclusions 

In  this  work  we  have  simulated  three  channel  flows  using  the  suboptimal-control  code 
developed  by  Nicoud  et  al.  (2001),  after  some  modifications  to  improve  its  performance. 
We  have  seen  that  the  suboptimal  control  code  is  not  able  to  reproduce  the  wake  compo- 
nent of  the  mean  velocity  profile  of  real  channels.  The  magnitude  of  the  error  depends  on 
the  parameter  a which  weights  the  energy  of  the  control  in  the  cost  function.  Decreasing 
this  parameter  decreases  the  error  in  the  wake,  but  degrades  the  agreement  of  the  velocity 
fluctuation  intensities  near  the  wall.  A simplified  analysis  of  two  numerical  experiments 
with  different  values  of  a suggests  that  the  problem  with  the  velocity  profile  near  the 
channel  centerline  is  intrinsic  to  the  simulation  procedure,  and  not  a consequence  of  the 
boundary  condition  algorithm.  This  result,  together  with  the  reorganization  of  the  flow 
structures  away  from  the  wall,  raises  the  question  of  whether  it  is  possible  to  control  the 
whole  flow  just  by  acting  on  the  wall. 

There  are  also  open  questions  about  the  effects  of  a on  the  mean  velocity  profile  and 
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Figure  8.  One-dimensional  premultiplied  velocities  spectra  at  y+  = 150,  normalized  in  wall 
units  and  plotted  as  functions  of  the  wavelengths  A*  (a,  c)  and  Az  (6,  d).  (a)  and  (fc),  streamwise 
component;  (c)  and  ( d)  spanwise  component. , DNS;  o , case  SI;  • , case  S2;  A , case  Rl. 

on  the  root  mean  squared  velocity  fluctuations,  where  the  influence  of  this  parameter 
seems  to  be  more  important.  More  work  is  needed  in  that  direction. 

We  have  shown  that  the  suboptimal  code  is  able  to  match  a synthetic  profile  for  a 
rough  turbulent  channel,  and  that  the  trend  of  the  changes  in  the  root-mean-squared 
velocity  fluctuation  profiles  agrees  with  the  expected  results.  It  is  not  clear  whether  this 
agreement  is  due  to  a change  in  the  physics  of  the  wall,  or  just  a secondary  effect  of  a lower 
mean  velocity  gradient.  However,  the  information  gathered  from  the  spectra  and  from 
the  p.d.f.’s  point  to  the  latter  explanation,  as  there  are  no  clear  differences  between  the 
structures  near  the  wall  in  the  smooth  and  the  rough  cases,  except  for  the  intensities.  The 
structure  of  the  wall  stresses  introduced  by  the  control  algorithm  bear  little  resemblance 
to  the  corresponding  physical  quantities. 

There  are  in  all  the  present  simulations  a spurious  organization  of  the  near-wall  span- 
wise  velocity  into  large  diagonal  modes,  which  is  most  likely  due  to  some  instability  of 
the  control  procedure,  but  which  also  disappears  a few  grid  points  away  from  the  wall. 
Its  analysis  also  requires  further  work. 
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